#Set options and load packages
options(stringsAsFactors=F)
rm(list=ls())
library(multiwayvcov)
library(lmtest)
library(lubridate)
library(boot)
library(zoo)
library(fixest)
library(PanelMatch)

#Read in data files
issuanceData<-read.csv("issuanceData.csv")
extensionData<-read.csv("extensionData.csv")
tradeVisaData<-read.csv("tradeVisaData.csv")

sink("log.txt", append = F, type = "output")

#Results Reported in the Main Text

#Produce Figure 1: Visa Issuance Before/After DMC Closure
#Isolate only those units that switch status
constantSwitchers<-issuanceData[issuanceData[,"Closure"]>=2008&issuanceData[,"Closure"]<=2016,]
#Limit to three years before/after switch in status
constantSwitchers<-constantSwitchers[abs(constantSwitchers[,"timeSince"])<=3,]
#Aggregate to quarterly level
qAgg<-aggregate(logIssuance~timeSince,FUN=mean,data=constantSwitchers)


#Plot the data
#Save to file
pdf(file="Figure1.pdf",width=7,height=5)
plot(qAgg,xlab="Years Before/After DMC Closure",ylab="Average Logged Visa Issuance",main="Visa Issuance before and after DMC Closure")
#Draw on 0 line and average lines
abline(v=0)
lines(c(-3,-0.25),c(mean(qAgg[qAgg[,"timeSince"]<0,"logIssuance"]),mean(qAgg[qAgg[,"timeSince"]<0,"logIssuance"])),lty=2)
lines(c(0.25,3),c(mean(qAgg[qAgg[,"timeSince"]>0,"logIssuance"]),mean(qAgg[qAgg[,"timeSince"]>0,"logIssuance"])),lty=2)
dev.off()
#Produce Table 3 (Regression of Logged Issuance by Country-Quarter)
#Model 1
feols(log(Issued+1)~closed | Nationality, cluster = ~Nationality+Quarter,data=issuanceData)
#Model 2
feols(log(Issued+1)~closed | Nationality + Quarter, cluster = ~Nationality+Quarter,data=issuanceData)
#Model 3
feols(log(Issued+1)~closed | Nationality^Year + Quarter, cluster = ~Nationality+Quarter,data=issuanceData)
#Model 4
feols(log(Issued+1)~closed | Nationality + Quarter^Region, cluster = ~Nationality+Quarter,data=issuanceData)
#Model 5
feols(log(Issued+1)~closed | Nationality^Year + Quarter^Region, cluster = ~Nationality+Quarter,data=issuanceData)


#Produce Table 4 (Regression of Logged Issuance by Country-Quarter Split by Visa Type)
#Model 1
feols(log(Issued.Visitor+1)~closed +onshoredV | Nationality + Quarter, cluster = ~Nationality+Quarter,data=issuanceData)
#Model 2
feols(log(Issued.Visitor+1)~closed +onshoredV | Nationality + Region^Quarter, cluster = ~Nationality+Quarter,data=issuanceData)
#Model 3
feols(log(Issued.Study+1)~closed +onshoredS | Nationality + Quarter, cluster = ~Nationality+Quarter,data=issuanceData)
#Model 4
feols(log(Issued.Study+1)~closed +onshoredS | Nationality + Region^Quarter, cluster = ~Nationality+Quarter,data=issuanceData)
#Model 5
feols(log(Issued.WF+1)~closed +onshoredWF | Nationality + Quarter, cluster = ~Nationality+Quarter,data=issuanceData)
#Model 6
feols(log(Issued.WF+1)~closed +onshoredWF | Nationality + Region^Quarter, cluster = ~Nationality+Quarter,data=issuanceData)


#Produce Table 5 (Regression of Logged Extension Grants by Country-Quarter)
#Model 1
feols(log(Granted.Extensions+1)~closed | Nationality,cluster=~Nationality + Quarter, data = extensionData)
#Model 2
feols(log(Granted.Extensions+1)~closed | Nationality + Quarter,cluster=~Nationality + Quarter, data = extensionData)
#Model 3
feols(log(Granted.Extensions+1)~closed | Nationality^Year + Quarter,cluster=~Nationality + Quarter, data = extensionData)
#Model 4
feols(log(Granted.Extensions+1)~closed | Nationality + Region^Quarter,cluster=~Nationality + Quarter, data = extensionData)
#Model 5
feols(log(Granted.Extensions+1)~closed | Nationality^Year + Region^Quarter,cluster=~Nationality + Quarter, data = extensionData)



##Models Reported in the Appendix##
##Note To Replicator: Some of these models, especially the models using the Sun-Abraham Correction in Tables 6 and 10 and are computationally intensive and will take a fair amount of time to run.  Note that Figure 2 (the event study plot) depends on running the Table 6 models first, but the rest of these models can be run independently##

#The sunab function expects to have always treated units with their treamnt period set to the first period in the data and wants periods as dates.  Reformat the data frame to conform to this expectation

issueSA<-issuanceData
issueSA[,"Closure"]<-as.yearqtr(issuanceData[,"Closure"],format="%Y %B")
issueSA[issueSA[,"Closure"]==1970,"Closure"]<-2005
issueSA[,"Quarter"]<-as.yearqtr(issueSA[,"Quarter"],format="%Y Q%q")


#Produce Table 6 (Regression of Logged Issuance by Country-Quarter Using Sun-Abraham Correction)
#Model 1
sa1<-feols(log(Issued+1)~sunab(Closure,Quarter) | Nationality, cluster = ~Nationality+Quarter,data=issueSA)
summary(sa1,agg='ATT')
#Model 2
sa2<-feols(log(Issued+1)~sunab(Closure,Quarter) | Nationality + Quarter, cluster = ~Nationality+Quarter,data=issueSA)
summary(sa2,agg='ATT')
#Model 4 (note there is no Model 3)
sa4<-feols(log(Issued+1)~sunab(Closure,Quarter) | Nationality + Quarter^Region, cluster = ~Nationality+Quarter,data=issueSA)
summary(sa4,agg='ATT')

#Produce Figure 2 (Event Study Plot)
pdf(file="Figure2.pdf",width=7,height=5)
iplot(sa1,xlim=c(-29,29),main="Estimated Effect of Closure on Logged Issuance")
dev.off()

#Produce Table 7 (Replication of Table 3 Using Issuance Rate)
#Model 1
feols(Issued/(Issued+Refused) ~ closed | Nationality, cluster = ~Nationality+Quarter,data=issuanceData,weights=sqrt(issuanceData[,"Issued"]+issuanceData[,"Refused"]))
#Model 2
feols(Issued/(Issued+Refused) ~ closed | Nationality + Quarter, cluster = ~Nationality+Quarter,data=issuanceData,weights=sqrt(issuanceData[,"Issued"]+issuanceData[,"Refused"]))
#Model 3
feols(Issued/(Issued+Refused) ~ closed | Nationality^Year + Quarter , cluster = ~Nationality+Quarter,data=issuanceData,weights=sqrt(issuanceData[,"Issued"]+issuanceData[,"Refused"]))
#Model 4
feols(Issued/(Issued+Refused) ~ closed | Nationality + Quarter^Region , cluster = ~Nationality+Quarter,data=issuanceData,weights=sqrt(issuanceData[,"Issued"]+issuanceData[,"Refused"]))
#Model 5
feols(Issued/(Issued+Refused) ~ closed | Nationality^Year + Quarter^Region , cluster = ~Nationality+Quarter,data=issuanceData,weights=sqrt(issuanceData[,"Issued"]+issuanceData[,"Refused"]))


#Produce Table 8 (Regression of Logged Applications by Country-Quarter)

#Model 1
feols(log(Applications+1)~closed | Nationality, cluster = ~Nationality+Quarter,data=issuanceData)
#Model 2
feols(log(Applications+1)~closed | Nationality + Quarter, cluster = ~Nationality+Quarter,data=issuanceData)
#Model 3
feols(log(Applications+1)~closed | Nationality^Year + Quarter, cluster = ~Nationality+Quarter,data=issuanceData)
#Model 4
feols(log(Applications+1)~closed | Nationality + Quarter^Region, cluster = ~Nationality+Quarter,data=issuanceData)
#Model 5
feols(log(Applications+1)~closed | Nationality^Year + Quarter^Region, cluster = ~Nationality+Quarter,data=issuanceData)


#Produce Table 9 (Matching Estimate of ATT of Closure)
#Standard errors by bootstrapping, set a seed to ensure replicability
set.seed(2024)
pm.results1<-PanelMatch(lag=4, lead =0:4, time.id="period",unit.id="unit",treatment="closed",data=tradeVisaData,qoi="att",refinement.method="none",outcome.var="logIssuance")
pe1<-PanelEstimate(sets=pm.results1,data=tradeVisaData,pooled=T,confidence.level=0.999,number.iterations=10000)
summary(pe1)

pm.results2<-PanelMatch(lag=4, lead =0:4, time.id="period",unit.id="unit",treatment="closed",data=tradeVisaData,qoi="att",refinement.method="mahalanobis",outcome.var="logIssuance",covs.formula= ~lag(logIssuance,4))
pe2<-PanelEstimate(sets=pm.results2,data=tradeVisaData,pooled=T,confidence.level=0.999,number.iterations=10000)
summary(pe2)

pm.results3<-PanelMatch(lag=4, lead =0:4, time.id="period",unit.id="unit",treatment="closed",data=tradeVisaData,qoi="att",refinement.method="mahalanobis",outcome.var="logIssuance",covs.formula= ~Exports+Imports)
pe3<-PanelEstimate(sets=pm.results3,data=tradeVisaData,pooled=T,confidence.level=0.99,number.iterations=10000)
summary(pe3)

pm.results4<-PanelMatch(lag=4, lead =0:4, time.id="period",unit.id="unit",treatment="closed",data=tradeVisaData,qoi="att",refinement.method="mahalanobis",outcome.var="logIssuance",covs.formula= ~Exports+Imports +lag(logIssuance,4))
pe4<-PanelEstimate(sets=pm.results4,data=tradeVisaData,pooled=T,confidence.level=0.999,number.iterations=10000)
summary(pe4)


#Produce Table 10 (Regression of Logged Issuance by Country-Quarter Split bu Visa Type Using Sun-Abraham Correction)

#Model 1
sa101<-feols(log(Issued.Visitor+1)~sunab(Closure,Quarter) +onshoredV | Nationality + Quarter, cluster = ~Nationality+Quarter,data=issueSA)
summary(sa101,agg="ATT")
#Model 2
sa102<-feols(log(Issued.Visitor+1)~sunab(Closure,Quarter) +onshoredV | Nationality + Region^Quarter, cluster = ~Nationality+Quarter,data=issueSA)
summary(sa102,agg="ATT")
#Model 3
sa103<-feols(log(Issued.Study+1)~sunab(Closure,Quarter) +onshoredS | Nationality + Quarter, cluster = ~Nationality+Quarter,data=issueSA)
summary(sa103,agg="ATT")
#Model 4
sa104<-feols(log(Issued.Study+1)~sunab(Closure,Quarter) +onshoredS | Nationality + Region^Quarter, cluster = ~Nationality+Quarter,data=issueSA)
summary(sa104,agg="ATT")
#Model 5
sa105<-feols(log(Issued.WF+1)~sunab(Closure,Quarter) +onshoredWF | Nationality + Quarter, cluster = ~Nationality+Quarter,data=issueSA)
summary(sa105,agg="ATT")
#Model 6
sa106<-feols(log(Issued.WF+1)~sunab(Closure,Quarter) +onshoredWF | Nationality + Region^Quarter, cluster = ~Nationality+Quarter,data=issueSA)
summary(sa106,agg="ATT")

#Produce Table 11 (Alternative FE Specifications of Table 5)
#Model 1
feols(log(Issued.Visitor+1)~closed +onshoredV | Nationality, cluster = ~Nationality+Quarter,data=issuanceData)
#Model 2
feols(log(Issued.Visitor+1)~closed +onshoredV | Nationality^Year + Quarter, cluster = ~Nationality+Quarter,data=issuanceData)
#Model 3
feols(log(Issued.Visitor+1)~closed +onshoredV | Nationality^Year + Region^Quarter, cluster = ~Nationality+Quarter,data=issuanceData)
#Model 4
feols(log(Issued.Study+1)~closed +onshoredS | Nationality, cluster = ~Nationality+Quarter,data=issuanceData)
#Model 5
feols(log(Issued.Study+1)~closed +onshoredS | Nationality^Year + Quarter, cluster = ~Nationality+Quarter,data=issuanceData)
#Model 6
feols(log(Issued.Study+1)~closed +onshoredS | Nationality^Year + Region^Quarter, cluster = ~Nationality+Quarter,data=issuanceData)
#Model 7
feols(log(Issued.WF+1)~closed +onshoredWF | Nationality, cluster = ~Nationality+Quarter,data=issuanceData)
#Model 8
feols(log(Issued.WF+1)~closed +onshoredWF | Nationality^Year + Quarter, cluster = ~Nationality+Quarter,data=issuanceData)
#Model 9
feols(log(Issued.WF+1)~closed +onshoredWF | Nationality^Year + Region^Quarter, cluster = ~Nationality+Quarter,data=issuanceData)


#Produce Table 12 (Use of Extensions as a Control Variable)
#Model 1
feols(log(Issued+1)~ closed + log(Granted.Extensions+1)| Nationality,cluster=~Nationality+Quarter,data=extensionData)
#Model 2
feols(log(Issued+1)~ closed + log(Granted.Extensions+1)| Nationality +Quarter,cluster=~Nationality+Quarter,data=extensionData)
#Model 3
feols(log(Issued+1)~ closed + log(Granted.Extensions+1)| Nationality^Year +Quarter,cluster=~Nationality+Quarter,data=extensionData)
#Model 4
feols(log(Issued+1)~ closed + log(Granted.Extensions+1)| Nationality +Quarter^Region,cluster=~Nationality+Quarter,data=extensionData)
#Model 5
feols(log(Issued+1)~ closed + log(Granted.Extensions+1)| Nationality^Year +Quarter^Region,cluster=~Nationality+Quarter,data=extensionData)


sink()